#include "stdio.h"
#include "stdlib.h"
#include <iostream>
#include <sstream>
#include "curand.h"
#include <curand_kernel.h>
#include "algorithm"
using namespace std;
#include "CStemCell.h"
#include "System.h"
#include "main.h"
#include <iomanip>
#include "random"

#define VARSIGMA 0.1

extern struct IMD _MD;
extern default_random_engine eng;
extern uniform_real_distribution<double> uniform_dist;
extern normal_distribution<double> normal_dist;

int global_id = 1;

__global__ void print(double * d_sys_cells, size_t n_col, size_t n_row,int max_len = 10)
{
    // 打印矩阵 
    for(int i =0; i < max_len; i ++ )
    {
        for(int j = 0; j < NROW; j ++)
        {
            printf("%.2f ",d_sys_cells[i*NROW + j]);
        }
        printf("\n");
    }
}

CSystem::CSystem(int N0)
{
    sys_u = 0.2;               // Initialized the microenvironement. In this study, we only consider the situation with unchange microenviroment.
    sys_NumCell = N0;          // Number of total cell numbers
    sys_NumCellProlif = 0;
    sys_NumCellResting = N0;

    sys_par.readPar(_MD.cellpar);
    cudaMalloc(&d_sys_par,sizeof(double) * NUMMUTPAR );

    double h_sys_par[NUMMUTPAR];
    for(int i = 0; i < NUMMUTPAR;i ++)
    {
        h_sys_par[i] = sys_par.getPar(i);
    }
    cudaMemcpy(d_sys_par,h_sys_par,NUMMUTPAR,cudaMemcpyHostToDevice);

    h_sys_cells = (double*)malloc(INITCELL *NROW * sizeof(double));
    // std::cout <<  INITCELL *NROW * sizeof(double)<< std::endl;
    cudaError_t status;
    status = cudaMalloc((void **)&d_sys_cells, INITCELL *NROW * sizeof(double));
    // 步骤3: 检查错误
    if (status != cudaSuccess) {
        // 处理错误
        std::cerr << "cudaMalloc failed!(in d_sys_cells)" << std::endl;
        exit(1);
    }

    
    status = cudaMalloc((void **)&d_N0,sizeof(double));

    if (status != cudaSuccess) {
        // 处理错误
        std::cerr << "cudaMalloc failed!(in d_N0)" << std::endl;
        exit(1);
    }
    status = cudaMalloc((void **)&d_N1,sizeof(double));

    if (status != cudaSuccess) {
        // 处理错误
        std::cerr << "cudaMalloc failed!(in d_N1)" << std::endl;
        exit(1);
    }
    // 分配一个cpu上的内存
    
    sys_varsigma = _MD.varsigma;
    sys_mutationrate = _MD.mutationrate;

    for(int i = 0; i < INITCELL; i++)
    {
        if(i < N0){
            for(int j = 0; j < NUMMUTPAR ;j ++){
                h_sys_cells[i*NROW + j] = sys_par.getPar(j);
            }
            h_sys_cells[i*NROW + 20] = 0.0; //处于静息态
            h_sys_cells[i*NROW + 19] = global_id++; //设定细胞编号
            h_sys_cells[i*NROW + 21] = 0.0; //设定细胞年龄
            h_sys_cells[i*NROW + 22] =  uniform_dist(eng);       // Initial state of the epigenetic state x1
            h_sys_cells[i*NROW + 23] =  uniform_dist(eng)*0.1;    // Initial state of the epigenetic state x2
        }else{
            for(int j =0; j < NROW ;j ++)
            {
                h_sys_cells[i*NROW +j] = -1.0;
            }
        }    
    }
    cudaMemcpy(d_sys_cells,h_sys_cells,INITCELL *NROW * sizeof(double),cudaMemcpyHostToDevice);
}

CSystem::~CSystem()
{
    cudaFree(d_sys_par);
    cudaFree(d_sys_cells);
    cudaFree(d_N0);
    cudaFree(d_N1);
    d_sys_par = nullptr;
    d_sys_cells = nullptr;
    d_N0 = nullptr;
    d_N1 = nullptr;
}

// 通过老细胞,得到新细胞
__device__ void  setNextCell(double * oldcell,double * newcell,int idx ,double mutationrate,int size = INITCELL)
{

    // 首先新细胞的id
    newcell[19] = 1 ;
    // 新细胞的profq
    newcell[20] = 0;
    // 新细胞的年龄
    newcell[21] = 0;
    // 为每个线程获取一个唯一的随机数生成器
    curandState localState;


    // // 初始化随机数生成器的状态
    curand_init(0, idx, 0, &localState);
    // 更新新细胞的表格遗传状态
    for(int j = 0; j < NUMVAR ;j++)
    {
        double phi;
        double a, b;
        double z;
        double x1 = oldcell[22],x2 = oldcell[23];
        switch(j)
        {
            case 0:
                phi = 0.08 + oldcell[13] * pow(oldcell[4] * x1, 1.8)/(1+pow(oldcell[4] * x1, 1.8));
                break;
            case 1:
                phi = 0.08 + oldcell[14] * pow(oldcell[5] * x2, 1.8)/(1+pow(oldcell[15] * x2, 1.8));
                break;
            default:
                phi = 0;
        }
        a = oldcell[8]  * phi;
        b = oldcell[8] * (1-phi);
        // 定义beta分布
        double x,y;

        // auto g_dist_1 = std::gamma_distribution<double>{a,1};
        // x = g_dist_1(eng);
        // auto g_dist_2 = std::gamma_distribution<double>{b,1};
        // y = g_dist_2(eng);
        // 更新新细胞
        newcell[22+j] = curand_uniform(&localState);
    }
    // 更新细胞参数
    for(int j = 0; j < NUMMUTPAR ;j ++)
    {
        newcell[j] = oldcell[j];
    }

    int pos =(int)(NUMMUTPAR * curand_uniform(&localState));
    double varofpar;
    if( curand_uniform(&localState) < mutationrate)  // Perform gene mutation with change in parameters.
    {
        varofpar = oldcell[pos];
        varofpar = varofpar * exp( VARSIGMA* curand_uniform(&localState) ) / exp((VARSIGMA * VARSIGMA) / 2);
        oldcell[pos] = varofpar;
    }
}


// 计算两个向量的距离
__device__ double vectorNorm(double *a, double *b, int size = NUMMUTPAR) {
    double sum = 0.0;
    for (int i = 0; i < size; ++i) {
        double diff = a[i] - b[i];
        sum += diff * diff;
    }
    return sqrt(sum);
}

__global__ void kernel(double * d_sys_cells, size_t n_col, size_t n_row,double u,int N0, double * sys_par,double dt,double mutationrate)
{
    int idx =  blockIdx.x * blockDim.x + threadIdx.x;
    if(idx < n_col && d_sys_cells[idx * NROW + 20] < 0)
    { 
        curandState localState;
        // 初始化随机数生成器的状态
        curand_init(0, idx, 0, &localState);
        //第一步,计算动力学参数.
        double beta0,beta;
        double theta;
        double x1 = d_sys_cells[idx*NROW + 22],x2 = d_sys_cells[idx*NROW + 23];
        beta0 = d_sys_cells[idx*NROW + 6];
        theta = d_sys_cells[idx*NROW + 12] * (1.0 + d_sys_cells[idx*NROW + 16] * pow(d_sys_cells[idx*NROW + 3] * x2, 6.0)/(1.0 + pow(d_sys_cells[idx*NROW + 3] * x2, 6.0)));
        beta = dt*beta0 * (1.0/(1.0 + N0/theta)) * ((d_sys_cells[idx*NROW + 0] * x1 + pow(d_sys_cells[idx*NROW + 1] * x1,6.0))/(1+pow(d_sys_cells[idx*NROW + 2] * x1,6.0))); 

        double kappa1,kappa2,dist,kappa;
        kappa1 = d_sys_cells[idx*NROW + 9] * 1.0/(1.0 + pow(d_sys_cells[idx*NROW + 7] * x1, 6.0));
        dist = vectorNorm(d_sys_cells,&d_sys_cells[idx * NROW]); //
        kappa2 = d_sys_cells[idx*NROW + 17] *(pow(dist,4))/(pow(dist,4)+pow(d_sys_cells[idx*NROW + 18],4));
        kappa = dt*(kappa1 + kappa2);

        double g,mu;
        g = 3.0 * pow(x2, u) * pow(1-x2, 1 - u);
        mu =dt*d_sys_cells[idx*NROW + 10]/(1 + d_sys_cells[idx*NROW + 15] * exp(g));
        double _rand;
        if ( d_sys_cells[idx*NROW + 20] == 0) {
            
            _rand = curand_uniform(&localState); //产生0-1之间的随机数
            if (_rand < kappa) { //分化,移除模拟
                d_sys_cells[idx*NROW + 20] = -1;
            }else if (_rand < kappa + beta )
            {
                d_sys_cells[idx*NROW + 20] = 1;
                d_sys_cells[idx*NROW + 21] = 0.0;
            }else {
            }
        }else if (d_sys_cells[idx*NROW + 20] == 1) {
            _rand = curand_uniform(&localState);
            if (_rand < mu) //细胞凋亡
            {
                d_sys_cells[idx*NROW + 20] = -1;
            }
            else {
                if (d_sys_cells[idx*NROW + 21] < d_sys_cells[idx*NROW + 11])//继续增殖
                {
                    d_sys_cells[idx*NROW + 21]  += dt;
                } else  //细胞分裂
                {
                    for(int number = 0; number < 2; number ++){
                        for(int k = 0; k < n_col;k++)
                        {
                            if(d_sys_cells[k*NROW + 20] < 0) {
                                int pos;
                                pos = k;
                                setNextCell(&d_sys_cells[idx*NROW],&d_sys_cells[pos*NROW],idx,mutationrate);
                            }
                            if(k == n_col)
                            {
                                printf("错误\n");
                            }
                        }
                    }
                    d_sys_cells[idx*NROW+20] = -1;
                }
            }
        }

    }
}
__global__ void upDateSys(double * d_sys_cells, size_t n_col, size_t n_row, double * N0 ,double * N1 )
{ // N0 表示环境中的细胞数量, N1 表示环境中处于静息态的细胞数量
   N0[0] = 0;
   N1[0] = 0;
   for(int i = 0; i < n_col ; i++)
   {
     if( fabs(d_sys_cells[i*NROW + 20] -0) <1e-5){
        N0[0] ++ ;
        N1[0] ++ ;
     }else if(fabs(d_sys_cells[i*NROW + 20] -1) <1e-5){
        N0[0] ++ ;
     }
   }
   //printf("%f ,%f",N0[0],N1[0]);
}

bool CSystem::SystemUpdate(double t)
{  
    //t时刻的对每一个细胞进行更新
     // 定义线程块大小和网格大小
     int blockSize = 2048; // 每个线程块 256 个线程
     int numBlocks = (INITCELL + blockSize - 1) / blockSize; // 计算需要的线程块数量
     
     // 调用核函数
     kernel<<<numBlocks, blockSize>>>(d_sys_cells,INITCELL,NROW,sys_u,sys_NumCellResting,d_sys_par,_MD.dt,sys_mutationrate);
    //开辟空间,用于记录环境中的数据
    
    // 更新环境中的变量.这里调用一个线程
    //print<<<1,1>>> (d_sys_cells,INITCELL,NROW);
    upDateSys<<<1,1>>>(d_sys_cells,INITCELL,NROW,d_N0,d_N1);
    cudaDeviceSynchronize(); // 等待GPU完成所有工作

    //cudaMemcpy(&sys_NumCell,d_N0,sizeof(double),cudaMemcpyDeviceToHost);

    //cudaMemcpy(&sys_NumCellResting,d_N1,sizeof(double),cudaMemcpyDeviceToHost);
    
    //std::cout << sys_NumCell << " 静息态:" << sys_NumCellResting << std::endl;
    
    if(sys_NumCell==0)
    {
        cout << " This text means cell_lists'size == 0\n";
        cout <<"当前时间: "<< t << endl;
        return(0);
    }

    if(sys_NumCell > MAXCELL ){
        std::cout << "细胞数量超过MAXCELL\n,终止计算";
        exit(0);
    }
    return(1);

}

void CSystem::readMutationrate(std::string filename)
{
    std::ifstream file(filename);
    std::string line;

    // 读取文件的每一行
    while (std::getline(file, line)) {
        // 使用stringstream来分割行
        std::istringstream iss(line);
        std::string key, value;
        if (std::getline(iss, key, '=')) {
            iss >> value;

            // 根据键填充结构体
            if (key == "varsigma") {
                this->sys_varsigma = std::stod(value);
            } else if (key == "mutationrate") {
                this->sys_mutationrate = std::stod(value);
            }
        }
    }
}

void CSystem::Run()
{
    double t;
    int step = 0;

    // string fn;
    // fn = _MD.mdcrd + string(".dat");

    // std::ofstream file(fn.c_str());
    // file << std::fixed << std::setprecision(2);
    // if (!file.is_open()) {
    //     std::cerr << "无法打开文件" << std::endl;
    //     exit(0);
    // }


    for (t=0; t<=_MD.T1; t=t+_MD.dt)
    {
        step=step+1;
        if(SystemUpdate(t))
        {
            //file << t << " " << sys_NumCell << std::endl;
            if((step-1)%_MD.ntpx==0){
                //OutPutSys(t);
                std::cout << "progress " << t/(_MD.T1) << std::endl;
            }
        }
        else{
            printf("Tumor cells cleaned at day %4.2f\n",t);
            break;
        }
    }
    //file.close();
}


void CSystem::OutPutSys(double t)
{
    // string filename;
    // int u = int(t);
    // filename =  _MD.mdcrd + string("-")  + std::to_string(u) + string(".dat");

    // std::ofstream file(filename.c_str());
    // if (!file.is_open()) {
    //     std::cerr << "无法打开文件" << std::endl;
    //     exit(0);
    // }
    // file << std::fixed << std::setprecision(6);
    // auto ite = d_sys_cells.begin();
    // for(; ite!=d_sys_cells.end(); ite++)
    // {
    //     // Edit to change you output items.
    //     file << ite->_cellid <<" " << ite->_X[0] <<" "<< ite->_X[1] << " " << ite->_ProfQ << " ";

    //     for (int i=0; i< NUMMUTPAR; i++)
    //     {
    //         if(i == 10 || i == 17 || i == 18){
    //             file << std::fixed << std::setprecision(6) << ite->_par.getPar(i)  <<" ";
    //         }else{
    //             file << std::fixed << std::setprecision(3) << ite->_par.getPar(i) << " " ;
    //         }
    //     }
    //     file << "\n";
    // }
    // file.close();
}
